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Abstract 

Finite  difference  schemes  for  the  evaluation  of  first  and  second  derivatives  are  presented. 
These  second  order  compact  schemes  were  designed  for  long-time  integration  of  evolution  equa¬ 
tions  by  solving  a  quadratic  constrained  minimization  problem.  The  quadratic  cost  function 
measures  the  global  truncation  error  while  taking  into  account  the  initial  data.  The  resulting 
schemes  are  applicable  for  integration  times  fourfold,  or  more,  longer  than  similar  previously 
studied  schemes.  A  similar  approach  was  used  to  obtain  improved  integration  schemes. 
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1  Introduction 


The  simulation  of  hyperbolic  partial  differential  equations  often  requires  long-time  integration. 
The  physical  phenomena  described  by  these  equations  typically  possess  a  range  of  space  and  time 
scales,  turbulent  fluid  flow  is  a  common  example.  Accurate  numerical  simulation  of  this  type  of 
processes  requires  proper  representation  of  all  the  relevant  physical  scales  in  the  numerical  model. 
These  requirements  lead  recently  to  new  interest  in  Fade  approximations  also  known  as  compact 
finite  difference  schemes  (7). 

Compact  finite  difference  schemes  had  long  been  known  and  used  in  numerical  analysis  [1,2,  3j. 
They  offer  a  means  of  obtaining  high  order  approximations  to  differential  operators  using  narrow 
stencils.  This  is  achieved  by  treating  the  sought  derivatives  as  unknowns  and  solving  a  system  of 
equations  for  them.  Typically,  the  resulting  matrices  are  tridiagonal  or  pentadiagonal,  which  can 
be  efficiently  solved  . 

In  [7]  a  class  of  highly  accurate  compact  schemes  for  first,  second  and  higher  derivatives  were 
presented  and  analyzed.  A  notion  of  resolving  efficiency  was  introduced  which  should  measure 
the  accuracy  with  which  the  finite  difference  approximation  represents  the  exact  solution  over 
the  full  range  of  length  scales  that  can  be  reali.’ed  on  a  given  grid.  This  criterion  was  then  used 
to  compare  various  schemes  and  motivated  the  design  of  a  new  class  of  schemes,  the  so-called 
schemes  with  spectral  like  resolution.  These  are  fourth  order  pentadiagonal  systems  with  seven 
points  stencils.  The'r  improved  resolution  characteristics  were  obtained  by  giving  up  on  high 
formal  accuracy;  instead,  requiring  that  the  symbol  of  the  discrete  difference  operator  should 
agree  with  the  differential  operator  at  three  prescribed  high  frequencies.  However,  the  resolving 
efficiency  is  a  too  crude  measure  as  it  aissumes  that  aU  frequencies  occur  with  similar  magnitude  in 
the  initial  data.  In  the  present  paper  it  is  shown  that  for  problems  with  various  initial  conditions 
these  schemes  are  far  from  optimal. 

A  different  measure  for  evaluating  finite  difference  schemes  is  the  L2  norm  of  the  local  trunca¬ 
tion  error.  This  measure  which  takes  into  account  the  Fourier  components  present  in  the  solution 
and  their  amplitude  was  applied  in  [9]  to  design  explicit  time  marching  schemes  (i.e.,  discretiz¬ 
ing  time  and  space  simultaneously)  by  solving  analytically  constrained  minimization  problems 
with  quadratic  cost.  This  error  measure  seems  more  adequate  for  comparing  difference  schemes. 
However,  the  simultaneous  treatment  of  time  and  space  results  in  very  complex  optimization 
problems.  The  generalization  of  the  approach  of  [9]  to  problems  in  higher  dimensions  requires 
solving  nonlinear  constrained  optimization  over  a  large  set  of  parameters.  It  is  also  hard  to  apply 
that  approach  to  compact  schemes.  This  complexity  and  the  use  of  analytic  rather  than  numerical 
methods  makes  the  suggested  approach  impractical. 

In  [4]  a  heuristic  derivation  was  done  by  minimizing  the  weighted  error  {in  the  Fourier  space) 
of  the  discrete  and  continuous  operators. 

The  present  paper  uses  the  same  cost  function  as  [9]  with  several  important  differences.  First, 
improved  bounds  on  the  truncation  error  were  derived.  These  enabled  us  to  treat  time  and 
space  discretizations  separately.  This  greatly  simplifies  the  minimization  problem  by  reducing 
it  to  two  lower  dimensional  problems.  Further  simplification  is  obtained  by  optimizing  each 
partial  derivative  separately,  instead  of  approximating  the  whole  differential  operator  as  was  done 
in  [9].  These  reductions  of  problem  complexity  resulted  in  a  simple  and  general  approach  to 
synthesis  of  discretization  schemes.  It  enabled  us  to  design  highly  accurate  compact  difference 
schemes  and  integration  formulas  for  various  operators  and  initial  data.  The  resulting  second 
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order  approximations  proved  to  be  robust  to  perturbations  in  the  spectrum  of  the  initial  data, 
exhibiting  resolution  superior  to  other  known  schemes. 

The  organization  of  the  paper  is  as  follows.  In  Section  2  Fourier  analysis  is  used  to  obtain 
bounds  on  the  truncation  error.  In  section  3  approximations  to  derivatives  are  presented,  for  the 
first  and  second  derivatives  and  first  derivative  at  mid  ceU  points.  Appendices  A-C  fist  coefficients 
for  these  derivatives  for  various  stencils  and  initial  conditions.  Improved  time  integration  schemes 
are  developed  in  section  4,  and  their  coefficients  are  listed  at  Appendix  D.  Section  5  discusses 
generalization  of  the  present  approach  to  more  complex  problems.  Numerical  results  are  presented 
in  section  6.  Concluding  remarks  are  made  in  section  7. 

2  Bounds  on  the  truncation  error 

The  application  of  Fourier  analysis  for  the  design  and  evaluation  of  finite  difference  schemes  can 
be  found  in  many  sources,  e.g.,  [9,  10].  In  [12]  the  use  of  Fourier  analysis  in  the  numerical 
approximation  of  hyperbolic  problems  is  extensively  discussed. 

In  the  following  section  bounds  on  the  L2  norm  of  the  error  in  the  discrete  solution  are  derived, 
accounting  for  the  effect  of  discretization  both  in  space  and  time.  These  estimates  are  used  in 
subsequent  sections  to  design  highly  accurate  schemes. 

Consider  a  linear  constant  coefficient  partial  differential  equation  with  periodic  boundary 
conditions  of  the  form  : 

I  =  i.  (2.1) 

u(i,0)  =  «o(a:)  (2.2) 

Further  assume  that  the  solution  of  equation  (2.2)  does  not  grow  in  time. 

The  discrete  analog  of  this  equation  can  be  written  as  : 

=  F(^A0<  (2.3) 

<  =  «o  (2.4) 

where  h  is  the  meshsize  in  space,  and  P(/i,At)  is  a  stable  finite  difference  approximation. 

We  would  like  to  bound  the  Lj  norm  of  the  error  in  the  discrete  solution,  for  the  initial  value 
uo,  given  by  : 

e\nAt-,  tto)  =  ]lti(nAt)  -  nlW^  =  -  P’*(/»,  Af)tioll*  (2.5) 

In  the  sequel  t  will  be  used  instead  of  nAt  to  simplify  notation.  The  Fourier  transform  of  eq. 
(2.5)  yields 

fi  2 

y  (2.6) 

^  _  P"(wh)|)*  ltio(a;h)]*du;  (2.7) 

Where  is  the  discrete  operator  approximating  L,  and  t,  are  their  corresponding  symbols. 
Thus,  the  space  and  time  discretizations  errors  can  be  bounded  separately. 
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Denote  L{ujh)  =  tf{{uh)  +  Li{u}h)  for  the  real  and  imaginary  parts  of  L{u!h),  respectively; 


and  use  a  similar  decomposition  for  L^{uh).  Then 

^L(uih)t  _  ^L^{wh.)t  _  ^L{wh)t  ^Ji  _  ^(L’Hu>h)-Li(u)k))t^{L)^(u/h)-LR{^h))t^  (2.8) 

The  assumption  that  the  solution  does  not  grow  in  time  implies  that 

<  I  t>0  (2.9) 

Assume  that  the  discrete  solution  does  not  grow  in  time,  either.  That  is 

^  j  >  Q  ^2.10) 

For  real  numbers  0,  a  with  q  <  0,  simple  geometric  considerations  yield  the  following  bound 

|1  _e'«e‘"|  <  |^j  +  |l -e“i  (2.11) 

Combining  bounds  (2.9)  and  (2.11)  and  assuming  L^{u!h)  -  Lft{u!h)  <  0  results  in  : 

<  |Lj(a;h)-i/(wh)|<  +  |l (2.12) 

<  \l'}{LJh)~  Li{uh)\t  +  \L'kiuh)- LR{ujh)\t  (2.13) 


If  Zfl(a;h)  —  Lfi(uh)  >  0,  this  bound  can  be  obtained  using  the  same  argument  when  factoring 

jjj  (2.8). 

Denote  by  €{t;  uq)  the  error  due  to  spatial  discretization  only,  when  the  initial  data  is  uq.  For 
a  final  time  T,  using  (2.13)  yields  : 

(|L)(a.h)  -  L}iwh)\  +  \L'k{<-^h)  -  LR{uh)\y  |uo(w/i)|'dw  (2.14) 

■f 

Therefore,  a  difference  scheme  minimizing  the  integral  (2.14)  with  respect  to  initial  value  uq 
will  better  resolve,  in  the  Lj  norm  sense,  the  frequencies  occurring  in  the  solution. 

The  time  integration  operator  satisfies  : 


'  j=i 

(2.15) 

Under  the  previous  assumptions 

<  1 

(2.16) 

lP(u;h)|  <  1 

(2.17) 

Therefore 

n— 1 

1  e^'‘<"'*)^'^‘P"“’--’(u;/i)l  <  Cn 

j=i 

(2.18) 

where  C  =  0(1).  Hence 

imixV 

(2.19) 
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Denote  by  e{t:uo)  the  error  due  to  time  discretization  only,  when  the  initial  value  is  uq-  For  a 
final  time  T,  the  bound  (2.19)  implies  : 

e\T,uo)  <  -  Piwh  f\uo(^h)\^d-^  (2.20) 

Combining  of  these  estimates  yields  a  bound  on  the  L2  norm  of  the  global  truncation  error  : 

e^(r;uo)  <  [|i^(u;fi)  -  L/(a;fi)l  +  |4(a;/i)  -  ifi(a;/i)|  +  -  P(a;fi)lj '  juo(‘^/i)|2du; 

(2.21) 

These  estimates  wiD  be  used  in  the  following  sections  to  design  improved  compact  schemes  and 
integration  formulas. 


3  Approximating  derivatives 

3.1  Approximation  of  the  first  derivative 

Consider  a  uniformly  spaced  mesh  whose  nodes  are  indexed  by  i  and  its  meshsize  is  given  by 
h  =  ^,  where  N  +  1  is  the  number  of  grid  points.  The  variable  at  node  i  is  i,  =  ih  and  the 
function  value  at  the  nodes,  /,  =  /(x,),  are  given  for  0  <  i  <  A.  An  approximation  /'  to  the  first 
derivative  ^(xi)  should  be  computed  as  a  linear  combination  of  the  function  values  at  neighboring 
grid  points.  Compact  finite  difference  schemes  regard  the  approximation  as  unknown  and  a 
system  of  equations  is  solved  to  approximate  the  first  derivative  at  all  nodes,  simultaneously. 
Thus,  unlike  in  finite  differences,  the  derivative  at  node  t  depends  on  function  values  at  aJl  other 
nodes. 

Following  [7]  we  use  approximations  of  the  the  form  : 


^fl-2  +  +  //  +  ^ 


/i+3  -  fi-3  ,  ,  fi+2 


6/1 


+  b- 


4h 


fi-2  ,  fi+\  -  fi-\ 
- f-  a- 


2h 


(3.1) 


A  second  order  approximation  can  be  obtained  by  adding  a  constrmnt  that  the  Taylor  expansion 
on  both  sides  should  agree  until  the  second  order  term.  The  following  relation  must  hold  : 


0"1-^"1"C  —  Id”  2oi  -j-  2/3 


(3.2) 


Higher  order  schemes  may  be  obtained  by  further  matching  the  next  terms  in  the  expansion  [7]. 
However,  in  this  paper  merely  second  order  accuracy  is  enforced. 

The  symbol  of  the  differentiation  operator  is  given  by  : 


Whereas  the  symbol  of  the  discrete  approximation  (3.1)  is 

.asin(u;/i)  +  I sin(2ti;/i)  +  I  sin(3<*;/i) 
^  ^  1  +  2acos(a;h)  +  2/3cos(2w/i) 


(3.3) 


(3.4) 
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In  view  of  the  bound  (2.14),  define  the  following  constrained  minimization  problem  which  its 
solution  should  yield  a  compact  scheme  with  improved  resolution  properties  : 

min  |Z^(wh)  -  L{u)h)\  |u(u;h)|^dui  (3.5) 

a,b,c,a,0  I  I 

under  the  constraint 

fl  +  6  +  c=  1  +  2q(  -f  2/3  (3.6) 

where  L{uih)  and  L^{u)h)  are  given  by  (3.3)  and  (3.4),  respectively. 

Although  the  problem  was  formulated  as  a  constrained  minimization  problem,  it  can  be  trans¬ 
formed  by  substitution  to  an  unconstrained  minimization  problem  on  a  reduced  set  of  parameters. 
Moreover,  setting  some  of  the  parameters  to  zero  further  reduces  the  dimension  of  the  problem. 
Since  tridiagonal  systems  of  equations  are  more  amenable  to  numerical  solution  than  pentadi- 
agonal  ones,  setting  /3  =  0  seems  a  plausible  choice.  Similar  considerations  might  suggest  using 
a  narrower  stencil  obtained  by  setting  c  =  0,  as  weU.  All  those  possibilities  are  presented  in 
Section  6,  and  several  sets  of  coefficients  for  different  initial  data  are  listed  in  Appendix  A. 


3.2  Approximation  of  the  second  derivative 

The  derivation  of  compact  schemes  for  the  second  derivative  proceeds  in  an  analogous  way  to  the 
first  derivative.  The  starting  point  is  an  approximation  of  the  form 


fi+3  -  2fi  -f  /i_3  ,  t  fi+2  -  2fi  +  fi-i  ,  _  fi+i  -  2fi  4-  fi-l 


9/i2 


+  b 


4h2 


+  0- 


(3.7) 

where  //'  is  the  approximation  to  the  second  derivative  at  node  i.  Matching  the  Taylor  series 
coefficients  on  both  sides  of  (3.7)  yields  condition  (3.2)  for  the  second  order  accuracy. 

The  symbol  of  the  second  derivative  is  given  by  : 


L{uh)  = 


(3.8) 


The  symbol  of  the  discrete  approximation  (3.7)  is  : 

fh/  L\  2a  (cos(u;h)  -  1)  +  ^  (cos(2u>/t)  -1)4-^  (cos(3u;/t)  —  1) 

^  ~  1  -I-  2a  cos(a;h)  +  2/3  cos(2wh) 

The  constrained  minimization  problem  which  solution  is  the  sought  scheme  can  be  formulated 

f 


as  : 


/n 

\t^{uh)  -  £,(a>/i)j  |u(u;/i)i’*dw 
■f  '  ' 


(3.9) 


(3.10) 


under  the  constraint 

o-f-b+c  =  1-f-  2q  +  2/3  (3.11) 

Now,  however,  L{u;h)  and  L^{uh)  are  given  by  (3.8)  and  (3.9),  respectively. 
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3.3  Approximating  first  derivative  on  a  cell  centered  mesh 

The  approximation  of  the  first  derivative  at  the  cell  centered  mesh  is  of  the  form  ; 


+  /i  +  OLfU\  +  /^/^2  -  ^ 


-  fi-^  fi+i  -  /,_3  /i+i  -  /,. 


(3.12) 


The  second  order  of  the  approximation  is  guaranteed  by  condition  (3.2). 
The  symbol  of  the  differentiation  operator  is  : 

L{uh)  =  iuh 

While  the  symbol  of  the  discrete  approximation  (3.12)  is  : 

=  ■2asin(if)  +  f  sin(3f^)  +  ^sin(^) 


1  +  2o!  cos(wft)  +  2/3  cos(2a>h) 


(3.13) 


(3.14) 


A  constrained  minimization  problem  of  the  same  type  as  in  the  previous  sections  was  formulated 
and  solved  for  these  symbols. 


4  Approximation  of  the  integration  operator 


The  design  of  integration  schemes  is  substantially  limited  by  the  stabibty  requirement  which 
renders  high  order  schemes  computationally  costly.  It  is  well  known  [5]  that  an  explicit  order 
Runge-Kutta  method,  with  k  >  4,  should  employ  at  least  A;+l  stages  or  function  evaluations.  This 
large  number  of  function  evzduations  makes  high  order  schemes  impractical.  Therefore,  efforts 
have  been  made  to  obtain  schemes  of  lower  order  with  improved  char  cteristics.  W’ithin  this 
approach,  the  free  variables  in  the  Runge-Kutta  schemes  were  set  to  yield  better  truncation  error 
[5]  or  extended  stability  region  [6].  A  generalization  of  this  idea  is  to  give  up  on  formal  accuracy 
in  order  to  obtain  better  approximation  of  the  wavenumbers  relevant  to  the  problem  solved. 

The  discrete  time  integration  of  linear  constant  coefficient  partial  differential  equation 


(4.1) 


amounts  to  approximation  of  the  exact  discrete  solution  e^^^uo-  Therefore,  the  integration  scheme 
may  be  written  as 

n 

Pn{L^At)='^ai{L^Aty  (4.2) 

«=o 

where  a,  may  depend  on  L^.  The  order  of  the  integration  scheme  is  determined  by  the  number 
of  first  terms  a;  which  agrees  with  the  Taylor  expansion  of  e^. 

The  derivation  of  the  integration  schemes  is  similar  to  that  of  derivative  discretization,  i.e., 
a  constrained  quadratic  optimization  problem  is  formulated  based  on  the  error  estimate  (2.20). 
The  solution  of  this  minimization  problem  yields  an  improved  integration  scheme.  However,  the 
derivation  of  integration  schemes  is  more  involved  than  the  generation  of  compact  schemes  since 
the  stability  condition  leads  to  a  nonlinearly  constrained  minimization  problem. 


Following  (2.20)  the  next  optimization  problem  is  defined  : 

'*  (4.3) 

.2 
n 

subject  to  the  constraints 

a«  =  t:  0  <  t  <  p  (4.4) 

t! 

\KiL>'{uh)At)\^  <  1 

where  is  the  discrete  approximation  of  L  and  p  is  the  order  of  the  n  stage  formula.  Condition 
(4.4)  can  be  treated  by  substitution,  but  the  stability  condition  requires  an  explicit  treatment. 

In  accordance  with  our  general  approach,  we  believe  that  second  order  formal  accuracy  suffices. 
It  remains  to  determine  the  number  of  stages  in  the  integration  formula.  This  should  be  chosen 
to  assure  that  the  error  in  space  and  time  discretizations  (2.14)  and  (2.20),  respectively,  will  be 
of  similar  magnitude.  In  the  present  work  five  stage  schemes  of  second  order  were  investigated, 
i.e.,  n  =  5  and  p  =  2.  Integration  formulas  were  obtained  for  optimized  seven  point  tridiagonal 
compact  schemes  approximating  the  first  derivative,  and  were  tested  for  the  advection  equation 
in  one  and  two  space  dimensions. 

An  important  feature  of  the  present  approach  is  that  once  a  feasible  minimum  has  been  found 
for  a  prescribed  initial  value  and  a  given  CFL  number,  the  generated  scheme  wiU  be  stable  for 
this  data.  This  might  enable  the  use  of  somewhat  larger  time  steps. 

5  Approximation  of  differential  operators 

The  method  introduced  in  the  previous  sections  for  generating  optimal  finite  difference  approxi¬ 
mations  for  derivatives  and  time  integration  schemes  for  linear  constant  coefficient  equations  can 
be  extended  to  more  general  equations.  These  ideais  can  be  easily  adapted  to  handle  with  similar 
efficiency  more  involved  problems. 

The  error  bounds  derived  in  Section  2  can  be  generalized  for  d  dimensional  problems;  noting 
that  the  same  proof  holds  for  the  d  dimensional  case  after  changing  the  integration  over 
to  multi  integration  over  the  box  This  suggests  that  approximation  of  the  differen¬ 

tial  equation  should  be  obtained  by  solving  constrained  optimization  problems  in  d  dimensional 
Fourier  space  for  a  large  set  of  parameters.  For  some  equations,  solving  Ihls  large  minimization 
problem  might  be  essential  to  achieve  accurate  schemes.  Quite  often,  though,  a  set  of  simpler 
minimization  problems  can  be  obtained  by  optimizing  each  partial  derivative  separately,  resulting 
in  highly  accurate  approximations. 

The  approach,  which  was  successfully  tested  in  the  present  paper,  divides  the  optimization 
process  into  two  stages.  In  the  first  stage  a  set  of  schemes  are  designed  for  a  large  enough 
variety  of  typical  initial  data  (e.g.,  Gaussians  with  different  parameters,  in  our  examples).  This 
precomputation  is  performed  once  and  its  results  are  used  in  subsequent  simulations.  In  the 
second  stage,  the  actual  simulation,  the  initial  data  uo  is  Fourier  transformed  to  obtain  uo. 
The  discretization  of  the  partial  derivatives  is  determined  by  approximating  uq  as  a  product  of 
one  dimensional  functions  for  which  optimized  schemes  werf  designed.  Each  partial  derivative 
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is  discretized  using  the  corresponding  one  dimensional  optimized  sciieme.  The  time  marching 
scheme  is  selected  from  the  set  of  scl  lie  corresponding  to  the  approximating  one  dimensional 
functions.  In  the  present  work  the  selection  was  done  by  computing  the  error  norm  of  each 
candidate  integration  scheme  when  applied  to  the  approximate  initial  data  with  the  already 
determined  discretizations,  and  selecting  the  minimum  norm  scheme.  I'his  computation,  too.  can 
be  done  prior  to  the  actual  simulation  for  a  large  set  of  approximated  initial  data,  riius.  the 
marching  scheme  selection  can  be  done  by  looking  up  in  a  precomputed  table.  The  robustness  of 
the  proposed  schemes  to  perturbations  in  the  initial  data  yields  this  optimization  very  efficient: 
as  can  be  seen  in  the  numerical  results  presented  in  Section  (i.  It  should  be  noted  that  tin- 
time  required  to  obtain  an  appropriate  scheme  using  this  approach  is  negligitiie  relative  to  the 
simulation  time. 

When  the  frequencies  present  in  the  solution  change  with  time.  e.g..  due  to  time  dependent 
source  term,  the  computation  of  the  optimizeil  schemes  should  be  repeated  once  a  large  ciinnilative 
change  ha,s  occurred.  Still,  the  relative  cost  of  of  this  computation  is  minimal. 

The  Fourier  transform  gives  the  energy  content  of  the  whole  initial  data.  It  may  occur,  that 
the  initial  data  is  smooth  at  some  regions  of  the  computational  domain  and  oscillatory  in  others, 
in  which  case  the  designed  approximation  will  give  good  performances  over  the  whole  domain. 
One  can  do  better  by  computing  a  different  scheme  for  each  region  and  using  a  weighted  sum 
of  the  resulting  schemes  near  region  boundaries.  This  requires  computing  the  Fourier  transform 
locally  in  each  region.  The  localization  to  a  particulai  region  can  be  achieve  by  multiplying  ito  by 
a  function  with  a  compact  support  which  encloses  the  region. 

In  some  cases,  systems  of  equations  may  be  treated  in  a  similar  way.  Look  first  at  a  one 
dimensional  first  order  system 

Ut  =  Aur  (5.6) 

u(i,0)  =  wo(i) 

where  z4  is  a  px  p  symmetric  matrix.  Let  A  =  AP  be  a  diagonal  matrix,  and  denote  v  =  Pu. 
The  discretization  of  the  system 


Vt  -  \Vx 

i;(2-,0)  =  Puo{x) 


(5.7) 


can  be  done  in  an  analogous  way  to  the  scalar  case,  except  for  the  time  marching  scheme  which 
should  be  chosen  from  a  set  of  candidate  schemes  (as  for  the  multidimensional  scalar  equations). 
Thus,  highly  accurate  discretization  of  the  system  (5.6)  can  be  achieved  by  first  discretizing  (5.7) 
and  using  the  identity  For  systems  in  higher  dimension 


Ut 


u{x,0)  = 


«=i 

Uo(x) 


du 

dx, 


(5.8) 


each  partial  derivative  should  be  optimized  separately.  In  this  ca,se,  we  require  that  all  .4/  will  be 
symmetric,  but  it  is  not  necessary  that  they  are  simultaneously  cliagonalizable. 

The  proposed  schemes  might  be  useful  for  nonlinear  equations,  as  well.  There,  one  should 
design  the  schemes  for  the  linearized  equation;  and  will  be  obliged  modify  them  once  a  large 
change  in  the  amplitude  of  the  wavenumbers  appearing  in  the  solution  occur. 
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6  Numeri*^  il  results 

6.1  Ap''  jximation  of  derivatives 

Thf  constrained  minimization  problem  (3.5)  for  the  space  discretization  can  be  easily  solved  by 
substitution  using  (3.2).  Differentiation  of  the  resulting  quadratic  form  provides  a  set  of  necessary 
conditions  holding  at  the  minimum.  This  nonlinear  system  can  be  solved  using  Newton  method, 
yielding  a  local  minimum.  Since  the  schemes  obtained  using  this  process  significantly  improve 
previously  known  schemes  [7],  no  attempts  were  made  to  find  the  other  zeroes  of  the  nonlinear 
system,  searching  for  better  minima. 

Three  types  of  schemes  were  studied  :  (a)  tridiagonal  with  five  points  stencil,  i.e.,  =  c  =  0  . 

(b)  tridiagonal  with  seven  points  stencil,  i.e.,  ii  —  0,  (c)  pentadiagonal  with  seven  points  stencil. 
The  initial  approximation  to  the  Newton  iteration  was,  typically,  a  compact  scheme  with  the  same 
structure,  taken  from  [7]. 

It  can  be  observed,  in  figures  I  and  5,  that  the  modules  of  symbol  of  the  optimized  penta¬ 
diagonal  scheme  for  the  first  and  second  derivative  is  larger  then  the  modulus  of  the  differential 
symbol.  This  error  is  exceedingly  larger  for  schemes  generated  to  appro.ximate  narrower  spectra. 
The  overshooting  occurs  in  the  highest  end  of  the  spectrum  for  wavenumbers  not  appearing  in 
the  solution.  However,  since  the  stability  of  a  scheme  is  determined  by  the  values  assumed  by 
[11],  this  type  of  scheme  is  applicable  only  with  small  CFL.  Moreover,  the  desired  ro¬ 
bustness  is  limited  by  this  phenomenon.  Therefore,  this  behavior  of  the  approximation  can  not 
be  ignored.  .4  possible  remedy  can  be  found  in  searching  for  other  minimizers  of  the  quadratic 
form.  Using  the  tridiagonal  scheme  as  initial  approximation  for  the  Newton  process  converged 
to  solut’ons  without  this  limiting  property  but  with  reduced  resolution,  similar  to  the  tridiago¬ 
nal  schemes.  Other  possible  directions,  e.g.,  further  looking  for  other  minima  or  penalizing  in 
the  cost  function  for  this  behavior  were  not  explored.  This  is  since  we  believe  that  for  practical 
applications  pentadiagonal  systems  are  too  costly  to  solve,  whereas  the  tridiagonal  schemes  offer 
similar  resolution  characteristics,  are  easier  to  solve  and  do  not  suffer  from  this  deficiency.  The 
pentadiagonal  scheme  are  given  mainly  for  theoretical  reasons  as  a  counterpart  to  the  spectral-like 
approximations. 

k  proper  appreciation  of  the  superiority  of  the  proposed  schemes  can  be  gained  by  integrating 
with  them  hyperbolic  equations  for  long  time,  provided  the  integration  process  introduces  only 
negligible  numerical  errors.  This  requirement  necessitates  either  using  high  order  integration 
schemes  or  employing  exact  integration,  as  was  done  in  the  present  work.  The  experiments 
described  in  next  subsections  clearly  demonstrate  the  superior  behavior  of  the  proposed  optimized 
schemes. 

6.1.1  Approximation  of  the  first  derivative 

Compact  finite  difference  schemes  were  designed  and  tested  for  initial  data  of  the  form  for 

several  values  of  q.  In  Figure  1,  the  symbols  of  schemes  corresponding  to  a  =  2  are  plotted,  as 
well  as  the  weighted  error 

|L{u;/i)- Z.'‘(a;/i)||ii(w/i)|  (6.9) 
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for  the  more  accurate  schemes.  The  coefficients  of  the  optimized  schemes  can  be  found  in  .4ppendix 
A.  The  coefficients  of  the  other  schemes  were  taken  from  [7],  For  scheme  (a)  the  coefficients  were 

a  =  =  0,a  =  y,6  =  =  0  (6.!0) 

The  coefficients  of  scheme  (c)  were 

The  coefficients  of  the  spectral-like  scheme  (e)  were 

a  =  0.5771439, >5  =  0.0896406, a  =  1.3025166,5  =  0.99355,c  =  0.03750245  (6.12) 


It  can  be  seen  that  each  optimized  scheme  better  approximates  the  differential  operator  than 
its  non-optirnized  counterpart.  In  Figure  1,  one  can  observe  that  although  the  symbol  of  the 
spectral-like  pentadiagonal  scheme  follows  the  differential  symbol  for  more  wavenumbers  than  the 
tridiagonal  scheme,  the  Lj  norm  of  truncation  error  of  tridiagonal  scheme  is  somewhat  smaller 
for  this  data.  This  can  be  explained  by  noting  that  the  error  of  the  tridiagonal  scheme  is  mainly 
in  the  high  frequencies  while  the  spectral-like  scheme  has  large  error  at  the  smoother  Fourier 
components  where  the  present  initial  data  has  more  energy.  The  spectral-like  scheme  attains 
better  resolution  at  the  expense  of  larger  error  in  lower  frequencies.  The  error  in  the  optimized 
schemes  is  significantly  smaller  than  in  their  counterparts.  More  precisely,  computing  the  error 
norms  reveals  that  the  error  in  the  tridiagonal  scheme  is  about  six  times  larger  than  in  the 
optimized  tridiagonal  scheme  while  error  norm  of  the  spectral-like  scheme  is  about  seventeen 
times  larger  than  in  the  optimized  pentadiagonal.  The  plot  of  the  absolute  value  of  the  error 
also  reveals  that  the  L2  norm  was  used  as  a  minimization  criteria.  This  can  be  seen  from  the 
several  sign  changes  of  the  error  of  the  optimized  schemes,  being  in  accordance  with  the  averaging 
property  of  the  choisen  norm. 

Figure  2  demonstrates  the  better  resolution  of  the  optimized  scheme  by  exact  integrating  in 
Fourier  space  on  a  32  points  grid  with  the  pentadiagonal  spectral-like  scheme  and  the  pentadiag¬ 


onal  optimized  scheme  the  equation 


dt  dx 


(6.13) 


cr  —  0.8  (a  being  the  CFL  number)  was  used.  It  is  shown  that  at  time  T  -  10000,  the  error  in  the 
solution  using  the  optimized  scheme  is  smaller  than  the  error  at  time  T  =  1000  when  using  the 
spectral-like  scheme.  This  suggests  that  the  optimized  scheme  can  be  used  for  integration  time 
at  least  ten  times  longer  than  the  spectral-like  scheme,  in  close  accordance  with  the  ratio  of  the 
error  norms. 

Figure  3  displays  the  scheme’s  robustness  to  perturbation  in  initial  condition.  The  solution 
integrated  with  the  optimized  scheme  far  better  approximates  the  exact  solution  than  the  one 
employing  the  spectral-like  approximation,  even  for  initial  data  different  from  the  one  it  was 
designed  to  resolve.  This  holds  for  both  smoother  and  more  oscillatory  initial  data.  Although 
those  examples  do  not  give  a  quantitative  view  on  the  relative  efficiency  of  the  schemes  for  those 
initial  data,  one  can  see  in  both  figures  that  by  the  time  the  solution  with  the  optimized  scheme 
developed  significant  error  the  error  in  the  one  corresponding  to  the  spectral-like  scheme  is  so 
large  it  no  longer  approximates  the  exact  solution. 
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Figurt'  4  shows  a  two  dimonsioiial  oquatiuii  which  dctiioiislralf’s  the  robustness  of  tlie  jjroposed 

schemes.  In  this  example  the  initial  data  %vas  taken  to  be  the  (iaussian  <“(-'*  + '>"4 1  rotated  at  an 

angle  of  Then  the  program  searched  for  initial  data  of  tiie  form  f  "('o +  1.  for  the  integers 

1  <  n\,ni  <  7  which  yielded  the  best  approximation  to  the  initial  data.  The  pentadiagonai 

*  2  2 

schemes  optimized  for  initial  data  and  were  then  used  to  compute  Uj-  and  iiy. 

respectively.  In  this  example  Ui  =  3  and  ~  ‘2.  The  resulting  semi  discrete  system  was  solved 
by  exact  integration  in  Fourier  space  on  a  32x32  grid.  The  plot  shows  a  cut  through  the  .solution  in 
the  X  direction  containing  the  maximum  point  of  the  solution.  While  the  solution  correspomling 
to  the  optimized  disretization  closely  approximates  the  exact  solution,  the  .solution  di.screlized 
with  the  spectral-like  scheme  bears  very  little  resemblance  to  the  exact  .solution. 


6.1.2  Approximation  of  the  second  derivative 

The  coefficients  of  compact  schemes  for  various  initial  conditions  of  the  form  c““*'  ran  be  found 
in  Appendix  B. 

Figure  5  plots  absolute  value  of  the  symbol  of  the  second  derivative  and  the  weighted  error, 
for  o  =  2.  The  parameters  of  the  optimized  schemes  can  be  foiind  in  Appendix  B.  The  coefficients 
of  the  other  schemes  were  taken  from  [7].  Scheme  (a)  is  given  by  : 

o  =  ^,,3  =  0,£[  = -^,6  =  ^,f  =  0  (6.14) 

The  coefficients  of  scheme  (c)  are  : 


a  = 


„  696-  1191a  ^ 

=  0,a  =  - rr: - -b 


2454a  -  294 


,  f  = 


38'  428  '  535 

The  coefficients  of  the  spectral-like  pentadiagonai  scheme  (e)  are  ; 


1179a  -  344 
2140 


(6.15) 


a  =  0.50209266,^3  =  0.05669169, n  =  0.21564935,6=  1.723322,  c  =  0.1765973  (6.16) 


It  can  be  seen  that  the  error  in  the  non  optimized  schemes  is  significantly  larger  than  in  the 
optimized  ones.  It  is  interesting  to  note  that,  again,  for  this  specific  data  the  L2  error  norm  of 
the  spectral-like  scheme  is  about  an  order  of  magnitude  larger  than  the  non  optimized  tridiagonal 
scheme.  This  phenomenon  suggests  that  the  re.solution  efficiency  is  a  poor  estimate  for  discretiza¬ 
tions  evaluation.  Computing  the  error  norms  reveals  that  the  error  in  the  optimized  tridiagonal 
scheme  is  about  seven  times  smaller  than  in  the  non  optimized  scheme,  whereas  the  error  in  the 
optimized  pentadiagonai  .scheme  is  seventy  times  smaller  than  the  spectral-Uke  scheme,  for  this 
given  data. 

The  efficiency  of  the  pentadiagonai  schemes  was  compared  by  integrating  the  wave  equation  : 


d'^u  d^U 
dt^  '■  dx^ 

This  equation  was  put  in  a  system  form  ; 


(6.17) 


(6.18) 
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This  system  was  exact  integrated  on  a  .'{2  points  grid  and  the  results  are  gi%eii  in  Figures  0  i 
demonstrating  the  improved  accuracy  of  the  optimized  scheme  and  its  robustness,  respectively. 
Figure  6  demonstrates  that  the  optimized  scheme  can  be  used  for  integration  time  37  times,  or 
more,  longer  than  the  spectral-like  scheme.  In  Figure  7  the  scheme  robustness  is  clearly  shown  for 
initial  data  smoother  or  more  oscillatory  than  the  data  for  which  the  scheme  was  design  ed.  In  both 
cases,  by  the  time  a  significant  error  occurs  in  the  solution  discretized  with  the  optimized  scheme 
the  solution  corresponding  to  the  spectral  like  scheme  totally  difTers  from  the  exact  solution. 

The  initial  solution  and  its  approximation,  for  the  two  dimensional  problem  in  Figure  8  were 
obtained  similarly  to  those  of  the  example  in  Figure  4.  While  the  solution  integrated  with  the 
optimized  scheme  closely  approximates  the  exact  solution,  it  is  hard  to  see  that  the  solution 
corresponding  to  the  spectral-like  scheme  indeed  approximates  the  same  problem. 


6.1.3  Mid  cell  approximation  of  the  first  derivative 

Appendix  C  lists  the  coefficients  of  schemes  designed  for  various  initial  data.  The  coefficients  of 
the  schemes  taken  from  [7]  are  listed  below.  Scheme  (a)  is  given  by  : 


a  =  =  0,a  =  5(3  -  2a),b=  ^(22a  -  l),c  =  0 

The  coefficients  of  scheme  (c)  are  : 

_  37950  -  397250  651 15a -3350  25669a -6114 

“~  354’^~  31368  ’  “  20912  “  62736 

The  coefficients  of  the  tenth  order  pentadiagonal  scheme  (f)  are  : 

96850  9675  683425  505175  69049 

“  ~  288529’^  "  577058’“  “  865587’^  “  577058’*^  "  11731174 


(6.19) 


(6.20) 


(6.21) 


The  standard  compact  schemes  give  very  good  resolution  in  this  form  (see  Figure  9),  thus, 
the  improvement  introduced  by  the  optimized  schemes  is  smaller.  Optimizing  the  Iridiagonal 
scheme  yields  a  6.5  smaller  error  norm  while  optimizing  the  pentadiagonal  scheme  yields  a  2.5 
times  smaller  norm.  In  this  case,  the  error  norm  of  the  optimized  tridiagona!  scheme  is  very  close 
to  that  of  the  non  optimized  pentadiagonal  scheme. 

An  interesting  option  suggested  by  this  approach  was  to  optimize  the  ^  operator  in  order  to 
get  the  best  approximation  for  for  given  initial  values.  This  has  been  done  for  the  tridiagonal 
scheme  which  was  u.sed  to  integrate  equation  (6.17).  It  was  compared,  in  Figure  10.  to  the 
tridiagonal  schem-  ;rom  [7]  where  both  are  used  to  approximate  the  second  derivative  in  solving 
the  one  dimensional  wave  equation  in  the  system  form  (6.18).  Again,  the  optimized  scheme  gives 
significantly  better  approximation. 


6.2  Approximate  time  integration 

The  constrained  minimization  (4.3)-(4.5)  was  solved  by  requiring  that  the  solution  will  touch 
the  stability  constraint  at  one  point  while  maintaining  global  stability  and  and  minimizing  the 
functional.  The  point  which  gives  best  result  was  found  by  exhaustive  search.  This  straightforward 
approach  yielded  the  local  minima  reported  in  this  paper.  .Somewhat  better  integration  schemes 
might  be  achieved  by  using  more  advanced  optimization  techniques  (8}. 
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According  to  the  general  approach  outlined  in  Section  5,  one  should  choose  an  integration 
scheme  which  yields  truncation  error  of  similar  magnitude  in  time  and  space.  Since  the  stability 
region  for  several  fifth  order  si.x  stages  explicit  Runge-Kutta  schemes  intersects  the  imaginary 
axis  only  in  a  smaU  neighborhood  of  the  origin  (5,  6]  disabling  time  marching  with  large  CFL,  the 
optimized  scheme  was  compared  with  the  four  stage  fourth  order  Runge-Kutta.  VVe  preferred  this 
five  stage  scheme  which  has  an  error  norm  about  five  times  larger  than  the  space  discretization  lo 
the  seventh  order  scheme  which  yields  an  error  norm  about  eleven  times  smaller  than  the  space 
discretization  because  of  its  lower  computational  cost. 

The  analysis  performed  in  Section  2  suggests  that  the  integration  operator  should  be  opti¬ 
mized  with  respect  to  the  spatial  discrete  operator  employed,  i.e..  to  minimize  j| /’( - 
{wa)A(|j^^  following  examples  is  the  tridiagonal  approximation  for  when  initial 

data  is  and  a  =  0.8.  .Appendix  D  contains  the  coefficients  of  integration  schemes  for  various 

initial  data  when  is  the  tridiagonal  scheme  optimized  for  the  same  initial  data  and  a  —  0.8. 

Figure  11  plots  the  real  and  imaginary  parts  of  versus  the  four  stage  fourth  order 

Runge-Kutta  and  the  improved  scheme.  The  norm  of  the  imaginary  part  of  the  error  was  reduced 
by  a  factor  of  31  while  its  real  part  was  reduced  by  merely  a  factor  of  2.3. 

Figures  12-13  shows  the  integration  of  the  advection  equation  with  those  scheme  on  a  32  points 
grid,  demonstrating  the  superior  efficiency  and  robustness  of  the  proposed  schemes.  In  Figure  12 
one  c<,n  see  that  the  optimized  scheme  can  be  used  for  at  least  three  times  longer  integration  time 
than  the  Runge-Kutta  scheme  applied  to  the  tridiagonal  scheme  from  [7].  The  computed  error 
norms  suggests  the  time  marching  error  is  dominant  in  all  examples. 

The  two  dimensional  example  in  Figure  14  summarizes  the  approach  suggested  in  this  work. 
It  compares  the  optimized  tridiagonal  scheme  combined  with  the  appropriate  integration  formula, 
to  fourth  order  Runge  Kutta  applied  to  non-optimized  tridiagonal  discretization.  Although  the 
analysis  in  Section  2  applies  only  to  constant  coefficient  problems,  this  example  shows  it  holds, 
heuristicly,  to  variable  coefficient  equations,  as  well.  The  initial  data  for  this  problem  was  obtained 
in  a  similar  manner  to  that  in  example  4.  However,  instead  of  comparing  the  solutions  computed 
on  the  32  x  32  grid  to  the  exact  solution,  they  are  compared  to  the  solution  on  a  64  x  64  grid 
which  was  integrated  with  the  optimized  scheme  designed  for  the  narrowest  computed  Gaussian 
(q  =  7).  The  initial  data  for  the  finer  grid  was  obtained  by  bibnear  interpolation  from  the  coarser 
grid.  It  can  be  seen  that  the  optimized  scheme  yields  significantly  more  accurate  solution. 

7  Conclusion 

A  simple  and  general  approach  for  the  design  of  finite  difference  approximation  of  derivatives 
and  integration  formulas  was  introduced.  It  was  used  to  design  compact  finite  difference  schemes 
for  derivatives  evaluation;  and  the  resulting  schemes  were  compared  to  previously  known  similar 
schemes.  The  guiding  line  has  been  to  improve  the  repre.sentation  of  the  range  of  wavenumbers 
appearing  in  the  physical  problem  being  .solved,  taking  into  account  their  relative  amplitudes. 
This  lead  to  an  measure  of  the  approximation.  The  resulting  schemes  combined  adaptivity  to 
the  specific  initial  data  by  the  nature  of  their  design  and  robustness  to  perturbations  in  the  initial 
data  .  The  improved  resolution  had  been  demonstrated  for  several  problems. 

A  similar  approach  was  tised  to  de.sign  improved  integration  .schemes,  taking  into  account  the 
spatial  discretization  as  well  as  the  initial  data. 
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The  approach  suggested  in  this  paper  for  optimizing  discrete  operators  can  be  similarly  applied 
to  higher, derivatives.  Its  applicability  to  more  general  and  complex  operators  should  be  further 
investigated. 

The  use  of  these  ideas  to  design  boundary  conditions  will  be  presented  elsewhere. 
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Coefficients  of  first  derivative  approximations  for  various  ini 
tial  conditions 


schemes  with  =  c  =  0 


Q  =  0.3793894912,  a  =  1.575573790,  b  = 


a  =  0.3534620453,  a  =  1.566965775,  b  = 


a  =  0.3461890571,  a  =  1.5633098070, b  = 


a  =  0.3427812069,  a  =  1.5614141543,  b  = 


Q  =  0.3408027739,  a  =  1.. 5602604992,  b  - 


a  =  0.3395099051,  a  =  1.5594855939,  b  - 


a  =  0.3385987444,  a  =  1.5589295176, b  = 


:  0  1832051925 


:  0.1399583152 


=  0.1290683071 


=  0.124148259 


=  0.121345048 


=  0.119534216 


=  0.1182679712 


a  =  0.4303030674,  a  = 


a  =  0.3991476265,0  = 


a  =  0.3904091387,0  = 


a  =  0.3863287472,  a  = 


a  =  0.3839604005,  o  = 


a  =  0.3824122042,0  = 


a  =  0.3813206436,0  = 


schemes  with  0  =  0 


:  1.5567577428, b  =  0.3451622238,  c  = 


:  1.5636386371, b  =  0.2563784492,  c  = 


:  1.5638887738,  b  =  0.2348222711,  c  = 


:  1.5637497712, b  =  0.2252138483,  c  = 


1.5635937780, 6  =  0.21976694619  ,  c 


1.5634617985, b  =  0.21625718276, c  = 


1.5633544597, 6  =  0.21380659696, c  = 


-0.0413138317 


-0.0217218334 


-0.0178927675 


-0.0163061252 


:  -0.0154399233 


:  -0.0148945794 


:  -0.0145197694 


tio 

general  schemes 

■1 

a  =  0.5779403671,  0  =  0.0890143475 
a  =  1.3030269541,  b  =  0.994883769, c  =  0.0359987066 

mu 

a  =  0.5801818925,  0  =  0.0877284887 
a  =  1.3058941939,  6  =  0.9975884963,  c  =  0.0323380724 

m 

o  =  0.5821143744,  0  =  0.0867224075 
a  =  1.3086733956,  b  =  0.9990906893,  c  =  0.0299094788 

■1 

a  =  0.5831688320,  0  =  0.0862000893 
a  =  1.3102698137,  b  =  0.9997174262,  c  =  0.0287506026 

^m 

tt  =  0.5838221871,  0  =  0.0858844217 
a  =  1.3112828763,  6  =  1.0000513827  ,  c  =  0.0280789585 

mm 

a  =  0.58426518608,  0  =  0.0856735831 
a  =  1.31197935750,  b  =  1.00025665126,  c  =  0.02764152958 

Q  =  0.58458494112,  0  =  0.08552292859 
o  =  1.31248665912,  b  =  1.00039487751,  c  =  0.02733420278 
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B  Coefficients  of  second  derivative  approximations  for  various 
initial  conditions 


schemes  with  =  c  =  0 


a  =  0.2285657609,  a  =  1.0139538409,  6  = 


a  =  0.2028150072,  a  =  1.05981351 70, b  = 


a  =  0.1952770765,  a  =  1.0716695072,  b  = 


a  =  0.1917151916,  a  =  1.0770076313, b  = 


a  =  0.1896428309, a  =  1.0800332355,  b  = 


a  =  0.18828772017,  a  =  1.0819792783,  b  ^ 
a  =  0.18733255632,  a  =  1.0833354275,  b  = 


0.4431776810 


0.3458164974 


0.3188846458 


0.3064227519 


0.29925242633 


0.29459616204 

0.29132968512 


Q  =  0.3125176074,  a: 


a  =  0.2702488609,  a  ■■ 


a  =  0.2580699154  ,  a 


a  =  0.2523894606,  a 


a  =  0.2491062584,  a  ■- 


a  =  0.2469677390,  a  ■- 


a  =  0.2454642305,  a  ■■ 


schemes  with 


0.7701351999,6  = 


0.8863525584, 6  = 


:  0.9170322739,6  = 


0.9308701065, 6  = 


0.9387256232, 6  = 


0.9437849227, 6  = 


0.9473144209, 6  = 


0  =  0 


0.9469577413,  c  = 


0.7065172637,  c : 


0.6425330979,  c 


0.6135153110,  r  : 


0.5969863585,  c  = 


0.5863166347,  c  = 


0.5788609571,  c: 


-0.0920577265 


-0.0523721002 


-0.043425.5409 


-0.0396064963 


-0.0374994649 


-0.0361660793 


-0.0352469171 


«0 

general  schemes 

mi 

a  =  0.5024750577, 0  =  0.0554440666 
a  =  0.2150536435,  b  =  1.7246523136,  c  =  0.1761322914 

mu 

a  =  0.5041582074,  0  =  0.0527585356 
a  =  0.2120465713,  b  =  1.7488409942,  c  =  0.15294.59205 

mm 

a  =  0.5053986368,  0  =  0.0512444502 
a  =  0.2112256102,6  =  1.7609579037,  c  =  0.1411026601 

mi 

a  =  0.5061009898,  0  =  0.0504756862 
a  =  0.2110263782,  6  =  1.7667867767,  c  =  0.1353401973 

mi 

a  =  0.5065435817,  0  =  0.0500170894 
a  =  0.2109783634,  6  =  1.7701652358  ,  c  =  0.1319777431 

■ 

a  =  0.5068465815,  0  =  0.0497133535 
a  =  0.2109761550,  6  =  1.7723629924,  c  =  0.1297807226 

mi 

a  =  0.5070666579  0  =  0.0494975852 
a  =  0.2109890794,  6  =  1.7739051293,  c  =  0.1282342776 
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C  Coefficients  of  mid  cell  approximation  of  the  first  derivative 
for  various  initial  conditions 


Mo 

.  schemes  with  /?  =  c  =  0 

WBM 

a  =  0.1824466564,  a  =  0.9847348088,  b  =  0.3801585039 

jBSM 

a  =  0.1621215357,0  =  1.002655871 1, 6  =  0.3215872003 

o  =  0.1560892225,  a  =  1.0076143702,  b  =  0.3045640747 

'JBSBM 

a  =  0.1532174394,  a  =  1.0099120.548,  b  =  0.2965228240 

a  =  0.1515399131,0  =  1.01 12348225,  6  =  0.2918450036 

^1 

a  =  0.1504402935,  a  =  1.0120939889,  b  =  0.288786-5980 

BBI 

«  =  0.1496639344,  a  =  1.0126967653,  b  =  0.2866311035 

Mo 

schemes  with  =  0 

a  =  0.2803531992,  a  =  0.8656018611,  b  0.7202754832,  c  =  -0.0251709460 

BSI 

a  =  0.2421691108,  a  =  0.9108711860,  b  =  0.5897758895,  c  =  -0.0163088538 

IBI 

a  =  0.2311768224,  a  =  0.9233491904,  b  =  0.5531540626,  c  =  -0.0141496081 

BBI 

Q  =  0.2260281312,  a  =  0.9290969691,  6  =  0.5361564844,  c  =  -0.0131971911 

BBi 

tt  =  0.2230456380,  a  =  0.9323967450,  b  =  0.5263571378  ,c  =  -0.0126626068 

a  =  0.2211004185,  a  =  0.9345368034,  b  =  0.5199847.302,  c  =  -0.0123206966 

a  =  0.2197316282,  a  =  0.9360368674,  6  =  0.51.55096846,  c  =  -0.0120832956 

Mo 

general  schemes 

a  -  0.3392424034, =  0.0126851467 
a  =  0.7880308119,  6  =  0.8956208871, c  =  0.0202034010 

m 

a  =  0.3364203680,  0  =  0.0159838314 
a  =  0.7894607720,  b  =  0.8790559502,  c  =  0.0362916767 

BH 

a  =  0.3359766282,  0  =  0.0164557610 
a  =  0.7895453413,  b  =  0.8768367139,  c  =  0.0384827231 

HH 

a  =  0.335834.5755, 13  =  0.0166014190 
aO.789.5561.5727,  6  =  0.87616875736,  c  =  0.03914707436 

u 

a  :=  0.33577201042,  (3  =  0.01666433833 
a  =  0.789.55722003,  6  =  0.87588406207  ,  c  =  0.03943141540 

gm 

a  =  0.33573907328,  0  -  0.01669706335 
a  =  0.78955658369,  b  =  0.87573722427,  c  =  0.03957846.531 

HH 

Q  =  0.33571963682,/?  =  1.67162152050 
a  =  0.78955572985,  b  =  0.87565178238,  c  =  0.03966419181 
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«0 

schemes  with  3  =  0,  designed  to  approximate 

e 

a  =  0.2949304593,  a  =  0.8473898079.  b  =  0.7718938474.  c  ^  -0.0294227367 

o  =  0.2482825125,  a  -  0.9037600128.  b  =  0.6104318128,  c  =  -0.0176268008 

r.  0.2349387889,  a  =  0.91908.59222.  b  ^  0., 56.56763 75 7.  r  =  -0.0148847202 

a  =  0.2287.3857.54,  a  ^  0.9260661646,  b  =  0..545H27700,  c  ^  -0.013701783s 

Kil 

Q  =  0.2251628777,  a  -  0.9300484196,  b  =  0.. 5333228985  .r  =  -0.013045.5627 

a  =  0.2228371850,  a  =  0.9326209622,  b  =  0..52.56822919,  r  =  -0.0126288841 

Q  =  0.2212037269,  a  =  0.9.344193325,  6  =  0..520329 10793,  c  =  -0.0123409866 

D  Coefficients  of  time  integration  scheme 


uo 

third  order  schemes  designed  for  a  —  0.9  having  gq  =  1 ,  Ci  =  1, 02  =  | 

Kil 

^3  =  0.166281  ,  =  0.0397196 , 05  =  0.0076705 

mm 

as  =  0.166407  ,  =  0.0409.525  ,  05  =  0.0074510 

BBI 

aa  =  0.1664488 , 04  =  0.04111513  ,  as  =  0.00739737 

EHI1 

as  =  0.1664805 , 04  =  0.04121264  ,  as  =  0.00736302 

Ki 

a3  =  0.1665028 , 04  =  0.04128218  ,  as  =  0.00733301 

EBI 

03  =  0.166.5207 , 04  =  0.04133150 , 05  =  0.00731074 
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Wavenumber— w  Wavenumber— w 


Figure  1:  Symbols  (left)  and  absolute  value  of  error  (right)  for  uo  =  (a)  Sixth  order 

tridiagonal  scheme  (/?  =  c  =  0)  (b)  Second  order  optimized  tridiagonal  scheme  {&  =  c  —  0)  (c) 
Eighth  order  tridiagonal  scheme  {(i  =  0)  (d)  Second  order  optimized  tridiagonal  scheme  (/?  =  0) 
(e)  Spectral-like  pentadiagonal  (f)  Optimized  pentadiagonal.  (g)  Exact  symbol.  Schemes  were 
optimized  for  uq  =  . 


T  =  1000.  T  =  10000. 


Fignr«  2;  Long  time  integration  of  the  equation  Uf  =  Uj,  uq  =  a  =  0.8.  (a)  Pentadiagonal 

scheme  optimized  for  uo  =  (b)  Spectral-like  pentadiagonal  scheme  (c)  Exact  solution. 
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T  »  1800.  T  =  15000. 


Figure  3:  Long  time  integration  of  the  equation  ut  —  <7  =  0.8.  Initial  solution  on  the  left 

figure  was  uo  =  ;  on  the  right  figure  it  was  Uq  —  e—*"  . 

(a)  Pentadiagonal  scheme  optimized  for  wo  =  (b)  Spectral-like  pentadiagonal  scheme  (c) 

Exact  solution. 
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T  =  6500. 


0  2  4  6 

X 


Figure  4:  Integration  of  the  equation  U(  =  +  Uy,  cr  =  0.8  using  pentadiagonal  schemes.  Initial 

solution  was  uq  —  rotated  at  an  angle  of  This  data  was  approximated  by  unrotated 

gaussian  (a)  Optimized  pentadiagonal  scheme  (b)  Spectral-like  pentadiagonal  scheme 

(c)  Exact  solution. 
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Wavenumber— ij 


Wavenumber— w 


Figure  5:  Symbol  (left)  and  absolute  value  of  error  (right)  for  ^  .  uq  =  (a)  Sixth  order 

tridiagonal  scheme  (/)  =  c  =  0)  (b)  Second  order  optimized  tridiagonal  scheme  (/3  =  c  =  0)  (c) 
Eighth  order  tridiagonal  scheme  {0  =  0)  (d)  Second  order  optimized  tridiagonal  scheme  {0  =  0) 
(e)  Spectral-like  pentadiagonal  (f)  Optimized  pentadiagonal.  (g)  Exact  symbol.  Schemes  were 
optimized  for  tlo  = 


T  =  200. 


T  =  7500. 


0.15 

0.10 

0.05 

0.00 

-0.05 

0  2  4  6 
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Figure  6;  Long  time  integration  for  -  it„,  «o  —  t.  cr  =  0.8.  (a)  Pentadiagonal  scheme 
optimized  for  uq  =  (b)  Spectral-like  pentadiagonal  scheme  (c)  Exact  solution. 
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T  =  3200. 


T  =  16000. 


Figure  7:  Long  time  integration  for  =  0.8.  Initial  solution  on  the  left  figure  wa.s 

uq  =  e""  ;  on  the  right  figure  it  was  tio  =  .  (a)  Pentadiagonal  scheme  optimized  for 

uo  =  (b)  Spectral-like  pentadiagonal  scheme  (c)  Exact  solution. 
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=  6500. 


Figure  8:  Long  time  integration  for  Uu  =  u^x  +  Uyy,  a  =  0.8  using  pentadiagonal  schemes. 
Initial  solution  on  the  left  figure  was  rio  =  rotated  at  an  angle  of  This  data 

was  approximated  by  unrotated  gaussian  (a)  Optimized  pentadiagonal  scheme  (b) 

Spectral-like  pentadiagonal  (c)  Exact  solution. 


Wavenumber— u 


Figure  9:  Symbol  for  mid-cell  discretizations  of  uo  =  (a)  Sixth  order  tridiagonal 

scheme  (/3  =  c  =  0)  (b)  Second  order  optimized  tridiagonal  scheme  (/?  =  c  =  0)  (c)  Eighth  order 
tridiagonal  scheme  (/3  =  0)  (d)  Second  order  optimized  tridiagonal  scheme  {3  =  0)  (e)  Tenth 
order  pentadiagonal  (f)  Optimized  pentadiagonal,  (g)  Exact  symbol.  Schemes  were  optimized  for 
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Figure  11:  Real  and  imaginary  parts  of  approximations  to  where  is  the  symbol 

of  the  tridiagonal  scheme  for  ^  optimized  for  uq  =  and  a  =  0.8.  (a)  Five  stage  scheme 

optimized  for  the  same  a  (b)  Fourth  order  Runge-Kutta  (c)  Exact  time  integration. 


T  -  ISO. 


T  -  450. 


Figure  12:  Integration  of  u*  =  Ux,  uo  =  er  =  0.8.  The  space  derivative  is  computed  using 

the  tridiagonal  compact  scheme  optimized  for  the  same  initial  date  and  a.  (a)  Five  stage  scheme 
optimized  for  this  scheme  and  CFL  (b)  Fourth  order  Runge  Kutta  (c)  Exact  time  integration 
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T  =  45. 


T  =  1000. 


Figure  13:  Integration  of  Uf  =  u^,  tr  =  0.8.  Left:  uo  =  Right:  uq  =  (a)  Five 

stage  scheme  optimized  for  this  scheme  and  CFL  (b)  Fourth  order  Runge  Kutta  (c)  Exact  time 
integration. 
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T  =  400. 


0  2  4  6 
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Figure  14:  Integration  of  ut  =  +  0.5(1  +  0.63tn{23ry))uy,  a  =  0.8  using  tridiagonal  schemes. 

Initial  solution  was  «o  =  rotated  at  an  angle  of  This  data  was  approximated 

by  unrotated  gaussian  (a)  Optimized  tridiagonal  scheme  and  optimized  marching 

scheme  (b)  Tridiagonal  scheme  integrated  with  fourth  order  Runge-Kutta.  (c)  A  fine  grid  solution 
(practically  exact) 
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